## define function for DE filtering: fig1-2, all supp fig
library(pheatmap)
library(RColorBrewer)
library(greekLetters)
library(tidyverse)
library(patchwork)
library(ggpubr)
library(ggpmisc)
source('DE_cells/scripts/DE_filtering_funcs.R')

# AID x 3 -------
aid_list <- list.files(path = "DE_cells/results", pattern = "ms|pss|^sle", full.names = TRUE) |>
  map(read_csv) |>
  map(\(x)select(x, -c(nCount_RNA:RNA_snn_res.0.8)))

# AID without soupx
aid_list <- list.files(path = "DE_cells/results/raw", pattern = "ms|pss|sle", full.names = TRUE) |>
  map(read_csv) |>
  map(\(x)select(x, -c(nCount_RNA:RNA_snn_res.0.8)))

## set patient & ctrl groups ----
aid_list |> map(dplyr::count, orig.ident)

### MS CSF & PBMC
aid_list[[1]] <- aid_list[[1]] |>
  mutate(group = str_remove(orig.ident, '[0-9]+') |>
           str_replace('PST|PTC', 'HC'))

aid_list[[1]] |> dplyr::count(group)

aid_list[[4]] <- aid_list[[1]] |>
  filter(str_detect(group, 'CSF'))

aid_list[[1]] <- aid_list[[1]] |>
  filter(str_detect(group, 'PBMC'))

### pSS
aid_list[[2]] <- aid_list[[2]] |>
  mutate(group = if_else(str_detect(orig.ident, 'HC'), 
                         'pSS_HC',
                         false = 'pSS'))

### SLE
aid_list[[3]] <- aid_list[[3]] |>
  mutate(group = if_else(str_detect(orig.ident, 'SLE'), 
                         'SLE',
                         false = 'SLE_HC'))

aid_list |>
  list_rbind() |>
  count(group, monaco_label) |>
  pivot_wider(names_from = group, values_from = n) |>
  write_csv('DE_cells/results/MS.table.4.aid.de.csv')

## DE cell frac for AID patient and ctrl -----
aid_de_cell <- aid_list |>
  map(calc_de_cell_frac, group, .progress = TRUE)|>
  map(mutate,
      disease = str_remove(group, '_HC'),
      disease = case_when(str_detect(disease, 'CSF|PBMC') ~ str_replace(disease, 'HC', 'MS'),
                          .default = str_glue('{disease}_PBMC')),
      group = case_when(str_detect(group, 'HC') ~ 'Healthy control',
                        .default = 'AID patient')) |>
  purrr::reduce(bind_rows) |> 
  mutate(type = case_when(str_detect(type, 'DE') ~ type,
                          .default = str_glue('{type}-')) |>
           fct_relevel("Ahmed DE", "Zxk DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |> 
  arrange(type) 

aid_pval <- aid_de_cell |>
  select(-fraction) |>
  mutate(group = case_match(group, 'AID patient' ~ 'AID', .default = 'HC'),
         total = total - count) |>
  pivot_wider(names_from = group, values_from = c(total, count)) |>
  rowwise() |>
  mutate(pval = t.test(c(rep(0,total_HC), rep(1,count_HC)),
                       c(rep(0,total_AID), rep(1,count_AID)))$p.value,
         psignif = case_when(pval < .001 ~ '***',
                             pval < .01 ~ '**',
                             pval < .05 ~ '*',
                             .default = ''))

aid_de_cell |> 
  mutate(type = fct_relevel(type, "Ahmed DE", "Zxk DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |>
  ggplot(aes(type, log10(fraction+1e-5), color = group, group = group)) + 
  geom_path() +
  geom_text(data = aid_pval, aes(label = psignif, x = type, y = -1.5), inherit.aes = FALSE) +
  theme_pubr() +
  labs_pubr() +
  expand_limits(y = c(-5, -1)) +
  scale_color_manual(values = c('red', 'blue')) +
  labs(x = 'DE cell defining criteria',
       y = 'log10(DE cell fraction + 1e-5)') +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5))+
  facet_wrap(~disease, scales = 'free')

aid_de_cell |>
  write_csv('DE_cells/results/aid_cell_cell_patient_HC.csv')

aid_de_cell <- read_csv('DE_cells/results/aid_cell_cell_patient_HC.csv')

## autoimmune DE per sample -----
aid_de_cell <- aid_list |>
  map(calc_de_cell_frac, orig.ident, .progress = TRUE)|>
  map2(c('MS_PBMC', 'pSS_PBMC', 'SLE_PBMC', 'MS_CSF'),
       \(x,y)mutate(x,
                    disease = y,
                    group = case_when(str_detect(orig.ident, 'MS|pSS|SLE') ~ 'AID patient',
                                      .default = 'HC'))) |>
  purrr::reduce(bind_rows) |> 
  mutate(type = case_when(str_detect(type, 'DE') ~ type,
                          .default = str_glue('{type}-')) |>
           fct_relevel("Ahmed DE", "Zxk DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |> 
  arrange(type) 

res <- aid_de_cell |>
  nest(data = -disease) |>
  mutate(plot = map2(data, disease, \(x,y){
    x |>
      mutate(type = fct_relevel(type, "Ahmed DE", "Zxk DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |>
      ggplot(aes(group, fraction, color = group)) +
      stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
      geom_jitter(height = 0, width = .1) +
      stat_compare_means(comparisons = list(c('AID patient','HC')), method = 't.test') +
      scale_y_continuous(expand = expansion(mult = c(0,.3))) +
      scale_color_manual(values = c('red','blue')) +
      theme_pubr() +
      facet_wrap(~type, scales = 'free') +
      ggtitle(y)}))

# pdf 8 * 10 inch
res$plot[[1]]
res$plot[[2]]
res$plot[[3]]
res$plot[[4]]

aid_de_cell |>
  write_csv('DE_cells/results/aid_HC_DE_per_sample.csv')

# covid x 2 ---------
cov_list <- list.files(path = "DE_cells/results/", pattern = "bgi|blish", full.names = TRUE) |>
  map(read_csv)

## set patient & ctrl groups ----
cov_list |> map(dplyr::count, orig.ident)

cov_list[[1]] <- cov_list[[1]] |>
  mutate(group = case_when(str_detect(orig.ident, 'COV') ~ 'Covid-19',
                           str_detect(orig.ident, 'Flu') ~ 'Flu',
                           .default = 'Healthy control'))

cov_list[[2]] <- cov_list[[2]] |>
  mutate(group = if_else(str_detect(orig.ident, 'covid'), 'Covid-19', 'Healthy control'))

cov_list |> map(dplyr::count, group)

cov_list[[1]] |>
  count(monaco_label, Stage) |>
  pivot_wider(names_from = Stage, values_from = n) |>
  write_csv('DE_cells/results/MS.table.1.bgi.de.csv')

cov_list[[2]] |>
  count(monaco_label, group) |>
  pivot_wider(names_from = group, values_from = n) |>
  write_csv('DE_cells/results/MS.table.2.blish.de.csv')

## DE cell frac for COVID patient and ctrl -----
cov_de_cell <- cov_list |>
  map(calc_de_cell_frac, group, .progress = TRUE)|>
  map2(c('BGI','Blish'),\(x,y)mutate(x, cohort = y)) |>
  purrr::reduce(bind_rows) |> 
  mutate(type = case_when(str_detect(type, 'DE') ~ type,
                          .default = str_glue('{type}-')) |>
           fct_relevel("Ahmed DE", "Zxk DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |> 
  arrange(type) 

cov_pval <- cov_de_cell |>
  select(-fraction) |>
  mutate(group = make.names(group),
         total = total - count) |>
  pivot_wider(names_from = group, values_from = c(total, count)) |>
  rowwise() |>
  mutate(pval = t.test(c(rep(0,total_Healthy.control), rep(1,count_Healthy.control)),
                       c(rep(0,total_Covid.19), rep(1,count_Covid.19)))$p.value,
         psignif = case_when(pval < .001 ~ '***',
                             pval < .01 ~ '**',
                             pval < .05 ~ '*',
                             .default = ''))

cov_de_cell |> 
  mutate(type = fct_relevel(type, "Ahmed DE", "Zxk DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |>
  ggplot(aes(type, log10(fraction+1e-5), color = group, group = group)) + 
  geom_path() +
  geom_text(data = cov_pval, aes(label = psignif, x = type, y = -1.5), inherit.aes = FALSE) +
  theme_pubr() +
  labs_pubr() +
  expand_limits(y = c(-5, -1)) +
  #scale_color_manual(values = c('red', 'blue')) +
  labs(x = 'DE cell defining criteria',
       y = 'log10(DE cell fraction + 1e-5)') +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5))+
  facet_wrap(~cohort, scales = 'free')

cov_de_cell |>
  write_csv('DE_cells/results/cov_cell_cell_patient_HC.csv')


## DE cell per sample in covid ----------
bgi_cov <- cov_list[[1]]

bgi_cov |>
  calc_de_cell_frac(orig.ident) |>
  mutate(group = case_when(str_detect(orig.ident, 'COV') ~ 'Covid-19',
                           str_detect(orig.ident, 'Flu') ~ 'Flu',
                           .default = 'HC')) |>
  ggplot(aes(group, fraction, color = group)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(comparisons = list(c('Covid-19','HC')), method = 't.test') +
  scale_y_continuous(expand = expansion(mult = c(0,.3))) +
  scale_color_manual(values = c('red','orange','blue')) +
  theme_pubr() +
  facet_wrap(~type, scales = 'free')

# per individual
bgi_cov |>
  mutate(individual = str_remove(orig.ident, '-D.+')) |>
  calc_de_cell_frac(individual) |>
  mutate(group = case_when(str_detect(individual, 'COV') ~ 'Covid-19',
                           str_detect(individual, 'Flu') ~ 'Flu',
                           .default = 'HC')) |>
  ggplot(aes(group, fraction, color = group)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(comparisons = list(c('Covid-19','HC')), method = 't.test') +
  scale_y_continuous(expand = expansion(mult = c(0,.3))) +
  scale_color_manual(values = c('red','orange','blue')) +
  theme_pubr() +
  facet_wrap(~type, scales = 'free')

# calc DE frac for cell types ------
aid_summary_celltype <- aid_list |>
  map(\(x)calc_de_cell_frac(x, monaco_label), .progress = TRUE) |>
  map2(c('MS_PBMC','pSS_PBMC','SLE_PBMC','MS_CSF'),
       \(x,y)add_column(x, disease = y)) |>
  purrr::reduce(bind_rows) |> 
  mutate(type = case_when(str_detect(type, 'DE') ~ type,
                          .default = str_glue('{type}-')) |>
           fct_relevel("Ahmed DE", "Zxk DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |> 
  arrange(type)

aid_summary_celltype <- aid_list[[4]] |>
  nest(.by = group) |>
  mutate(data = map(data, calc_de_cell_frac, monaco_label)) |>
  unnest(data) |>
  mutate(type = case_when(str_detect(type, 'DE') ~ type,
                          .default = str_glue('{type}-')) |>
           fct_relevel("Ahmed DE", "Zxk DE", 'CD79A-', 'CD79B-', 'CD79-', 'CD247-', 'CD3D-', 'CD3E-', 'CD3G-', 'CD3-')) |> 
  arrange(type)

aid_summary_celltype |>
  ggplot(aes(y = monaco_label, x = type, fill = fraction)) +
  geom_raster() +
  scale_fill_gradient(low = 'blue', high = 'red') +
  geom_text(aes(label = paste0(round(fraction*100, 1),'%')), color = 'yellow', size = 3) +
  theme_pubr() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
  labs(x = 'DE cell defining criteria',
       y = 'Cell types') +
  facet_wrap(~group)

## focus on MS_CSF -----
ms_csf <- aid_list[[4]] |>
  annotate_complex()

ms_csf |>
  filter(str_detect(monaco_label, 'Switched')) |>
  ggplot(aes(x = group, fill = (Zxk == TRUE))) +
  geom_bar() +
  theme_pubr() +
  scale_fill_manual(values = c('grey','red'), label = c('Non-DE cell', 'DE cell')) +
  labs(fill = 'Cell type',
       title = 'TCR+CD3+BCR+CD79+ criteria')

p1 <- last_plot()

ms_csf.cd3 <- aid_list[[4]] |>
  annotate_complex() |>
  leave_one_out('CD3') |>
  annotate_complex2()

ms_csf.cd3 |>
  filter(str_detect(monaco_label, 'Switched')) |>
  ggplot(aes(x = group, fill = (Zxk == TRUE))) +
  geom_bar() +
  theme_pubr() +
  scale_fill_manual(values = c('grey','red'), label = c('Non-DE cell', 'DE cell')) +
  labs(fill = 'Cell type',
       title = 'TCR+BCR+CD79+ criteria')

p2 <- last_plot()

p1 + p2 + plot_layout(guides = 'collect')+
  plot_annotation(theme = theme(legend.position = 'top'),
                  title = 'Switched memory B cells in CSF')

## examine Ig usage in DE vs non-DE switched memory B --------------
ms_csf.cd3 |>
  filter(str_detect(monaco_label, 'Switched') & group == 'MS_CSF') |>
  select(c(contains('IGH'), Zxk)) |>
  mutate(Zxk = case_when(Zxk == 0 ~ 'Non-DE cell', .default = 'DE cell')) |>
  pivot_longer(contains('IGH'), names_to = 'IGH genes', values_to = 'Expression') |>
  ggplot(aes(x = `IGH genes`, y = Expression, fill = `IGH genes`)) +
  geom_boxplot() +
  theme_pubr() +
  facet_wrap(~Zxk, nrow = 2)

ms_csf.cd3 |>
  filter(str_detect(monaco_label, 'B cells')) |>
  ggplot(aes(UMAP_1, UMAP_2, color = (Zxk == TRUE))) +
  geom_point()

aid_de_cell |>
  purrr::reduce(bind_rows) |>
  filter(type == 'Ahmed DE' | type == 'Zxk DE') |>
  select(-total) |>
  pivot_longer(c(count, fraction), names_to = 'Monaco_label', values_to = 'value') |>
  pivot_wider(names_from = 'group', values_from = value) |>
  unite('Monoca_label', 1:2) |>
  bind_rows(cell_type_count_t5) |>
  write_csv('DE_cells/results/table5_aids_monaco_labels.csv')

## find these cells' barcodes! ------
barcodes_ms_csf_de <- aid_list[[4]] |>
  fill_expr_mat() |>
  annotate_complex() |>
  select(.cell, Zxk)

barcodes_ms_csf_cd3_de <- aid_list[[4]] |>
  fill_expr_mat() |>
  annotate_complex() |>
  leave_one_out('CD3') |>
  annotate_complex2() |>
  select(.cell, Zxk)

count(barcodes_ms_csf_de, Zxk)

count(barcodes_ms_csf_cd3_de, Zxk)

write_csv(barcodes_ms_csf_de, 'DE_cells/results/barcodes_MS-CSF-DE.csv')
write_csv(barcodes_ms_csf_cd3_de, 'DE_cells/results/barcodes_MS-CSF-CD3-DE.csv')

## try visualize enrichment -----
de_cell_lou_tidy <- aid_de_cell |>
  purrr::reduce(bind_rows) |>
  mutate(disease_sample = str_extract(group, 'CSF|PBMC|pSS|SLE'),
         disease_sample = if_else(str_detect(group, 'CSF|PBMC'),
                                  str_glue('MS_{disease_sample}'),
                                  disease_sample),
         group = if_else(str_detect(group, 'HC'), 'control', 'disease'))

de_cell_lou_tidy |>
  mutate(pseudo_frac = (count + 0.1)/total) |>
  filter(disease_sample == 'MS_PBMC') |>
  ggplot(aes(type, count, fill = group)) +
  geom_col(position = 'dodge')

de_cell_aid_enrich <- de_cell_lou_tidy |>
  mutate(pseudo_frac = (count + .1)/total) |>
  pivot_wider(names_from = group, values_from = pseudo_frac, values_fill = 0) |>
  group_by(type, disease_sample) |>
  reframe(total_DE_count = sum(count),
          control = sum(control),
          disease = sum(disease),
          AID_DE_enrichment = disease / control)

de_cell_aid_enrich |>
  ggplot(aes(AID_DE_enrichment, type, fill = total_DE_count)) +
  geom_col() +
  facet_wrap(~disease_sample, scales = 'free')

## table 5
cell_type_count_t5 <- aid_list |>
  map(\(x) {
    x |>
      group_by(group) |>
      count(monaco_label) |>
      pivot_wider(names_from = group, values_from = n)
  }) |>
  purrr::reduce(full_join)

cell_type_count_t5

# plot zxk DE frac with 95 Conf.int -------
aid_de_cell <- read_csv('DE_cells/results/aid_cell_cell_patient_HC.csv')
cov_de_cell <- read_csv('DE_cells/results/cov_cell_cell_patient_HC.csv')

aid_de_conf <- aid_de_cell |>
  filter(type == 'Zxk DE') |>
  select(-fraction) |>
  mutate(total = total - count,
         group = case_match(group, 'Healthy control' ~ 'HC', 'AID patient' ~ 'AID')) |>
  pivot_wider(names_from = group, values_from = c(total, count)) |>
  rowwise() |>
  mutate(pval = t.test(c(rep(0,total_HC), rep(1,count_HC)),
                       c(rep(0,total_AID), rep(1,count_AID)))$p.value ,
         HC_es = t.test(c(rep(0,total_HC), rep(1,count_HC))) |> broom::tidy() |> select(estimate, conf.high),
         AID_es = t.test(c(rep(0,total_AID), rep(1,count_AID))) |> broom::tidy() |> select(estimate, conf.high)) |>
  select(2,7:9) |>
  pivot_longer(3:4) |>
  mutate(frac = value$estimate,
         conf.high = value$conf.high,
         name = case_match(name, 'AID_es' ~ 'patient', .default = 'control'), .keep = 'unused')

aid_conf_pval <- aid_de_conf |>
  filter(name == 'patient') |>
  mutate(group1 = 'patient',
         group2 = 'control',
         y.position = conf.high * 1.1,
         pval = signif(pval, digits = 3))

# save as 6*6 inch pdf
aid_de_conf |>
  ggplot(aes(name, frac, ymax = conf.high, ymin = frac, color = name)) +
  geom_col(fill = 'white', linewidth = 2) +
  geom_errorbar(width = .3) +
  facet_wrap(~disease, scales = 'free') +
  theme_pubr() +
  stat_pvalue_manual(aid_conf_pval, label = 'pval') +
  scale_y_continuous(expand = expansion(mult = c(0.1,0.1))) +
  labs(x = 'group', y = 'DE cell fraction', color = 'group') +
  scale_color_manual(values = c('blue','red'))

## covid-19 DE conf int plot -----
cov_de_conf <- cov_de_cell |>
  filter(type == 'Zxk DE') |>
  select(-fraction) |>
  mutate(total = total - count,
         group = case_match(group, 'Healthy control' ~ 'HC', 'Covid-19' ~ 'Cov', 'Flu'~'Flu')) |>
  pivot_wider(names_from = group, values_from = c(total, count), values_fill = 0) |>
  rowwise() |>
  mutate(pval = t.test(c(rep(0,total_HC), rep(1,count_HC)),
                       c(rep(0,total_Cov), rep(1,count_Cov)))$p.value ,
         HC_es = t.test(c(rep(0,total_HC), rep(1,count_HC))) |> broom::tidy() |> select(estimate, conf.high),
         Cov_es = t.test(c(rep(0,total_Cov), rep(1,count_Cov))) |> broom::tidy() |> select(estimate, conf.high)) |>
  select(!matches('count|total')) |>
  pivot_longer(contains('es')) |>
  mutate(frac = value$estimate,
         conf.high = value$conf.high,
         name = case_match(name, 'Cov_es' ~ 'Covid-19', 'Flu_es' ~ 'Flu', .default = 'control'), .keep = 'unused')

cov_conf_pval <- cov_de_conf |>
  summarise(conf.high = max(conf.high), .by = c(cohort,pval)) |>
  mutate(group1 = 'Covid-19',
         group2 = 'control',
         y.position = conf.high * 1.1,
         frac = y.position,
         pval = signif(pval, digits = 3))

# save as 6*6 inch pdf
cov_de_conf |>
  ggplot(aes(name, frac, ymax = conf.high, ymin = frac, color = name)) +
  geom_col(fill = 'white', linewidth = 2) +
  geom_errorbar(width = .3) +
  facet_wrap(~cohort, scales = 'free') +
  theme_pubr() +
  stat_pvalue_manual(cov_conf_pval, label = 'pval') +
  scale_y_continuous(expand = expansion(mult = c(0.1,0.1))) +
  labs(x = 'group', y = 'DE cell fraction', color = 'group',
       title = 'DE cells in Covid-19 patient PBMC') +
  scale_color_manual(values = c('blue','red'))

g1 <- last_plot()

## CRC DE conf int plot -----
chen2021 <- read_csv('DE_cells/results/chen2021-crc_de_meta.csv.gz')

crc_de_conf <- chen2021 |>
  filter(type == 'Zxk DE') |>
  select(-fraction) |>
  mutate(total = total - count) |>
  pivot_wider(names_from = tissue, values_from = c(total, count), values_fill = 0) |>
  rowwise() |>
  mutate(pval = t.test(c(rep(0,total_Normal), rep(1,count_Normal)),
                       c(rep(0,total_Primary), rep(1,count_Primary)))$p.value ,
         Normal_es = t.test(c(rep(0,total_Normal), rep(1,count_Normal))) |> broom::tidy() |> select(estimate, conf.high),
         Primary_es = t.test(c(rep(0,total_Primary), rep(1,count_Primary))) |> broom::tidy() |> select(estimate, conf.high)) |>
  select(!matches('count|total')) |>
  pivot_longer(contains('es')) |>
  mutate(frac = value$estimate,
         conf.high = value$conf.high,
         name = case_match(name, 'Primary_es' ~ 'tumor', .default = 'normal'), .keep = 'unused')

crc_conf_pval <- crc_de_conf |>
  summarise(conf.high = max(conf.high), .by = c(pval)) |>
  mutate(group1 = 'tumor',
         group2 = 'normal',
         y.position = conf.high * 1.1,
         frac = y.position,
         pval = signif(pval, digits = 3))

# save as 6*6 inch pdf
crc_de_conf |>
  ggplot(aes(name, frac, ymax = conf.high, ymin = frac, color = name)) +
  geom_col(fill = 'white', linewidth = 2) +
  geom_errorbar(width = .3) +
  theme_pubr() +
  stat_pvalue_manual(crc_conf_pval, label = 'pval') +
  scale_y_continuous(expand = expansion(mult = c(0.1,0.1))) +
  labs(x = 'group', y = 'DE cell fraction', color = 'group',
       title = 'DE cells in colon cancer tissue') +
  scale_color_manual(values = c('blue','red'))

g2 <- last_plot()

g1 + g2 + plot_layout(widths = c(1.5,1))

# perez SLE data --------
perez <- read_csv("DE_cells/results/perez_sle_de_meta.csv.gz")

perez |> dplyr::count(Status)

perez <- perez |>
  mutate(individual = str_c(ind_cov, '_', Status))

res_perez <- perez |>
  mutate(limit_de = CD247 & CD3E & CD3D & CD3G & CD79B & CD79A) |>
  group_by(Status) |>
  reframe(
    count = sum(limit_de),
    fraction = sum(limit_de) / n(),
    total = n()
  )

calc_conf <- function(target, other){
  seq <- c(rep(1,target), rep(0,other))
  t.test(seq) |>
    broom::tidy() |>
    select(estimate, conf.low, conf.high)
}

res_perez <- res_perez |>
  rowwise() |>
  mutate(conf = calc_conf(count, total-count)) |>
  unnest(conf)
  
res_perez |>
  mutate(Status = fct_relevel(Status, 'Healthy','Managed','Flare','Treated')) |>
  ggplot(aes(x = Status, y = fraction, ymin = conf.low, ymax = conf.high, color = Status)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .5) +
  scale_color_manual(values = c('blue','orange','red','purple')) +
  theme_pubr() +
  labs_pubr()

## by individual -----

ind_perez <- perez |>
  mutate(limit_de = CD247 & CD3E & CD3D & CD3G & CD79B & CD79A) |>
  group_by(individual) |>
  reframe(
    count = sum(limit_de),
    fraction = sum(limit_de) / n(),
    total = n()
  )

ind_perez <- perez |>
  dplyr::count(individual, Status, Sex, pop_cov, SLE_status, Age) |>
  select(-n) |>
  left_join(ind_perez)

ind_perez |>
  ggplot(aes(Status, fraction, color = Status)) +
  geom_jitter(width = .1, height = 0) +
  stat_summary(geom = 'crossbar')

ind_perez |>
  ggplot(aes(SLE_status, fraction, color = SLE_status)) +
  geom_jitter(width = .1, height = 0) +
  stat_summary(geom = 'crossbar') +
  stat_compare_means(comparisons = list(c('Healthy','SLE')), method = 't.test')

ind_perez |>
  ggplot(aes(Sex, fraction, color = Sex)) +
  geom_jitter(width = .1, height = 0) +
  stat_summary(geom = 'crossbar') +
  stat_compare_means(method = 't.test')

ind_perez |>
  ggplot(aes(pop_cov, fraction, color = pop_cov)) +
  geom_jitter(width = .1, height = 0) +
  stat_summary(geom = 'crossbar') +
  stat_compare_means()

ind_perez |>
  ggplot(aes(Age, fraction)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  theme_pubr() +
  labs_pubr() +
  stat_cor() +
  labs(y = 'DE cell fraction in PBMC')


res_perez |>
  ggplot(aes(type, log10(fraction+1e-5), color = group, group = group)) + 
  geom_path() +
  #geom_text(data = aid_pval, aes(label = psignif, x = type, y = -1.5), inherit.aes = FALSE) +
  theme_pubr() +
  labs_pubr() +
  expand_limits(y = c(-5, -1)) +
  #scale_color_manual(values = c('red', 'blue')) +
  labs(x = 'DE cell defining criteria',
       y = 'log10(DE cell fraction + 1e-5)') +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5))

res_perez |>
  filter(type == 'Zxk DE') |>
  ggplot(aes(group, fraction, fill = group)) +
  geom_col() +
  theme_pubr() +
  labs_pubr() +
  ggtitle('Fraction of IGH+IGL+CD79+TCRab/gd+CD3+ DE cells')

## per indv -----
indv_perez <- perez |>
  calc_de_cell_frac(ind_cov)

indv_perez <- perez |>
  dplyr::count(ind_cov,Status) |>
  right_join(indv_perez) |>
  filter(type == 'Zxk DE')

indv_perez |>
  ggplot(aes(Status, fraction, color = Status)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  theme_pubr() +
  labs_pubr() +
  ggtitle('Fraction of IGH+IGL+CD79+TCRab/gd+CD3+ DE cells')
